/*************************************************************************
*  2012-2020 Václav Šmilauer                                             *
*  2020-2025 Janek Kozicki                                               *
*                                                                        *
*  This program is free software; it is licensed under the terms of the  *
*  GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/

#pragma once
#include "common.hpp"
// classes dealing with exposing actual types, with many switches inside depending on template arguments

// helper function conv::val for converting strings to int, RealHP or ComplexHP
namespace conv {
template <typename Scalar, int Level = -1> inline typename boost::enable_if_c<std::is_same<int, Scalar>::value, int>::type val(const std::string& s)
{
	return boost::lexical_cast<int>(s);
}
template <typename Rr, int Level = ::yade::math::levelOfRealHP<Rr>>
inline typename boost::enable_if_c<::yade::math::isRealHP<Rr>, Rr>::type val(const std::string& s)
{
	return ::yade::math::fromStringRealHP<Rr>(s);
}
template <typename Cc, int Level = ::yade::math::levelOfComplexHP<Cc>>
inline typename boost::enable_if_c<::yade::math::isComplexHP<Cc>, Cc>::type val(const std::string& s)
{
	return ::yade::math::fromStringComplexHP<Cc>(s);
}
} // conv

// methods common for vectors and matrices
template <typename MatrixBaseT> class MatrixBaseVisitor : public py::def_visitor<MatrixBaseVisitor<MatrixBaseT>> {
	typedef typename MatrixBaseT::Scalar     Scalar;     // could be complex number
	typedef typename MatrixBaseT::RealScalar RealScalar; // this is the "real" (math) scalar
public:
	template <class PyClass> void visit(PyClass& cl) const
	{
		cl.def(py::init<MatrixBaseT>(py::arg("other")))
		        .def("__neg__", &MatrixBaseVisitor::__neg__<MatrixBaseT>)
		        .def("__add__", &MatrixBaseVisitor::__add__)
		        .def("__iadd__", &MatrixBaseVisitor::__iadd__)
		        .def("__sub__", &MatrixBaseVisitor::__sub__)
		        .def("__isub__", &MatrixBaseVisitor::__isub__)
		        .def("__eq__", &MatrixBaseVisitor::__eq__)
		        .def("__ne__", &MatrixBaseVisitor::__ne__)
		        .def("__mul__", &MatrixBaseVisitor::__mul__scalar<long>)
		        .def("__imul__", &MatrixBaseVisitor::__imul__scalar<long>)
		        .def("__rmul__", &MatrixBaseVisitor::__rmul__scalar<long>)
		        .def("isApprox",
		             &MatrixBaseVisitor::isApprox,
		             (py::arg("other"), py::arg("prec") = Eigen::NumTraits<Scalar>::dummy_precision()),
		             "Approximate comparison with precision *prec*.")
		        .def("rows", &MatrixBaseT::rows, "Number of rows.")
		        .def("cols", &MatrixBaseT::cols, "Number of columns.");
		visit_if_float<Scalar, PyClass>(cl);
		visit_fixed_or_dynamic<MatrixBaseT, PyClass>(cl);

		// reductions
		cl.def("sum", &MatrixBaseT::sum, "Sum of all elements.")
		        .def("prod", &MatrixBaseT::prod, "Product of all elements.")
		        .def("mean", &MatrixBaseT::mean, "Mean value over all elements.")
		        .def("maxAbsCoeff", &MatrixBaseVisitor::maxAbsCoeff, "Maximum absolute value over all elements.");
		// reductions only meaningful for non-complex matrices (e.g. maxCoeff, minCoeff)
		visit_reductions_noncomplex<Scalar, PyClass>(cl);
	};

private:
	template <class PyClass> static std::string name(PyClass& cl) { return py::extract<std::string>(cl.attr("__name__"))(); }
	// for dynamic matrices/vectors
	template <typename MatrixBaseT2, class PyClass, typename boost::enable_if_c<MatrixBaseT2::RowsAtCompileTime == Eigen::Dynamic, int>::type = 0>
	static void visit_fixed_or_dynamic(PyClass&)
	{
		//std::cerr<<"MatrixBaseVisitor: dynamic MatrixBaseT for "<<name(cl)<<std::endl;
		;
	}
	// for static matrices/vectors
	template <typename MatrixBaseT2, class PyClass, typename boost::disable_if_c<MatrixBaseT2::RowsAtCompileTime == Eigen::Dynamic, int>::type = 0>
	static void visit_fixed_or_dynamic(PyClass& cl)
	{
		//std::cerr<<"MatrixBaseVisitor: fixed MatrixBaseT for "<<name(cl)<<std::endl;
		cl.add_static_property("Ones", &MatrixBaseVisitor::Ones)
		        .add_static_property("Zero", &MatrixBaseVisitor::Zero)
		        .def("Random", &MatrixBaseVisitor::Random, "Return an object where all elements are randomly set to values between 0 and 1.")
		        .staticmethod("Random")
		        .add_static_property("Identity", &MatrixBaseVisitor::Identity);
	}
	template <typename Scalar, class PyClass, typename boost::enable_if<boost::is_integral<Scalar>, int>::type = 0> static void visit_if_float(PyClass&)
	{ /* do nothing */
	}
	template <typename Scalar, class PyClass, typename boost::disable_if<boost::is_integral<Scalar>, int>::type = 0> static void visit_if_float(PyClass& cl)
	{
		// operations with other scalars (Scalar is the floating type, long is the python integer type)
		// __trudiv__ is for py3k
		cl.def("__mul__", &MatrixBaseVisitor::__mul__scalar<Scalar>)
		        .def("__rmul__", &MatrixBaseVisitor::__rmul__scalar<Scalar>)
		        .def("__imul__", &MatrixBaseVisitor::__imul__scalar<Scalar>)
		        .def("__div__", &MatrixBaseVisitor::__div__scalar<long>)
		        .def("__truediv__", &MatrixBaseVisitor::__div__scalar<long>)
		        .def("__idiv__", &MatrixBaseVisitor::__idiv__scalar<long>)
		        .def("__itruediv__", &MatrixBaseVisitor::__div__scalar<long>)
		        .def("__div__", &MatrixBaseVisitor::__div__scalar<Scalar>)
		        .def("__truediv__", &MatrixBaseVisitor::__div__scalar<Scalar>)
		        .def("__idiv__", &MatrixBaseVisitor::__idiv__scalar<Scalar>)
		        .def("__itruediv__", &MatrixBaseVisitor::__idiv__scalar<Scalar>)
		        //
		        .def("norm", &MatrixBaseT::norm, "Euclidean norm.")
		        .def("__abs__", &MatrixBaseT::norm)
		        .def("squaredNorm", &MatrixBaseT::squaredNorm, "Square of the Euclidean norm.")
		        .def("normalize", &MatrixBaseT::normalize, "Normalize this object in-place.")
		        .def("normalized", &MatrixBaseT::normalized, "Return normalized copy of this object")
		        .def("pruned",
		             &MatrixBaseVisitor::pruned,
		             py::arg("absTol") = 1e-6,
		             "Zero all elements which are greater than *absTol*. Negative zeros are not pruned.");
	}
	// for fixed-size matrices/vectors only
	static RealScalar  maxAbsCoeff(const MatrixBaseT& m) { return m.array().abs().maxCoeff(); }
	static MatrixBaseT Ones() { return MatrixBaseT::Ones(); }
	static MatrixBaseT Zero() { return MatrixBaseT::Zero(); }
	static MatrixBaseT Random() { return MatrixBaseT::Random(); }
	static MatrixBaseT Identity() { return MatrixBaseT::Identity(); }

	static bool __eq__(const MatrixBaseT& a, const MatrixBaseT& b)
	{
		if (a.rows() != b.rows() || a.cols() != b.cols()) return false;
		return a.cwiseEqual(b).all();
	}
	static bool isApprox(const MatrixBaseT& a, const MatrixBaseT& b, const RealScalar& prec) { return a.isApprox(b, prec); }
	static bool __ne__(const MatrixBaseT& a, const MatrixBaseT& b) { return !__eq__(a, b); }

	// older compilers don't have std::void_t
	template <typename... Ts> struct make_void {
		typedef void type;
	};
	template <typename... Ts> using void_type = typename make_void<Ts...>::type;
	template <typename T, class = void> struct HasNegation : std::false_type {
	};
	template <typename T> struct HasNegation<T, void_type<decltype(std::declval<T>().operator-())>> : std::true_type {
	};
	// in older versions of Eigen AlignedVector3 has no negation operator. Following detection if it's present means that it will be used if a newer version of Eigen is used.
	template <typename MatrixBaseT2, typename boost::enable_if<HasNegation<MatrixBaseT2>, int>::type = 0> static MatrixBaseT __neg__(const MatrixBaseT& a)
	{
		return -a;
	};
	template <typename MatrixBaseT2, typename boost::disable_if<HasNegation<MatrixBaseT2>, int>::type = 0> static MatrixBaseT __neg__(const MatrixBaseT& a)
	{
		return Eigen::AlignedVector3<Scalar>(-a.x(), -a.y(), -a.z());
	};
	static MatrixBaseT __add__(const MatrixBaseT& a, const MatrixBaseT& b) { return a + b; }
	static MatrixBaseT __sub__(const MatrixBaseT& a, const MatrixBaseT& b) { return a - b; }
	static MatrixBaseT __iadd__(MatrixBaseT& a, const MatrixBaseT& b)
	{
		a += b;
		return a;
	};
	static MatrixBaseT __isub__(MatrixBaseT& a, const MatrixBaseT& b)
	{
		a -= b;
		return a;
	};

	template <typename Scalar2, typename boost::enable_if<std::is_convertible<Scalar2, const Scalar&>, int>::type = 0>
	static MatrixBaseT __mul__scalar(const MatrixBaseT& a, const Scalar2& scalar)
	{
		return a * scalar;
	}
	template <typename Scalar2, typename boost::disable_if<std::is_convertible<Scalar2, const Scalar&>, int>::type = 0>
	static MatrixBaseT __mul__scalar(const MatrixBaseT& a, const Scalar2& scalar)
	{
		return a * static_cast<Scalar>(scalar);
	}
	template <typename Scalar2, typename boost::enable_if<std::is_convertible<Scalar2, const Scalar&>, int>::type = 0>
	static MatrixBaseT __imul__scalar(MatrixBaseT& a, const Scalar2& scalar)
	{
		a *= scalar;
		return a;
	}
	template <typename Scalar2, typename boost::disable_if<std::is_convertible<Scalar2, const Scalar&>, int>::type = 0>
	static MatrixBaseT __imul__scalar(MatrixBaseT& a, const Scalar2& scalar)
	{
		a *= static_cast<Scalar>(scalar);
		return a;
	}
	template <typename Scalar2, typename boost::enable_if<std::is_convertible<Scalar2, const Scalar&>, int>::type = 0>
	static MatrixBaseT __rmul__scalar(const MatrixBaseT& a, const Scalar2& scalar)
	{
		return a * scalar;
	}
	template <typename Scalar2, typename boost::disable_if<std::is_convertible<Scalar2, const Scalar&>, int>::type = 0>
	static MatrixBaseT __rmul__scalar(const MatrixBaseT& a, const Scalar2& scalar)
	{
		return a * static_cast<Scalar>(scalar);
	}
	template <typename Scalar2, typename boost::enable_if<std::is_convertible<Scalar2, const Scalar&>, int>::type = 0>
	static MatrixBaseT __div__scalar(const MatrixBaseT& a, const Scalar2& scalar)
	{
		return a / scalar;
	}
	template <typename Scalar2, typename boost::disable_if<std::is_convertible<Scalar2, const Scalar&>, int>::type = 0>
	static MatrixBaseT __div__scalar(const MatrixBaseT& a, const Scalar2& scalar)
	{
		return a / static_cast<Scalar>(scalar);
	}
	template <typename Scalar2, typename boost::enable_if<std::is_convertible<Scalar2, const Scalar&>, int>::type = 0>
	static MatrixBaseT __idiv__scalar(MatrixBaseT& a, const Scalar2& scalar)
	{
		a /= scalar;
		return a;
	}
	template <typename Scalar2, typename boost::disable_if<std::is_convertible<Scalar2, const Scalar&>, int>::type = 0>
	static MatrixBaseT __idiv__scalar(MatrixBaseT& a, const Scalar2& scalar)
	{
		a /= static_cast<Scalar>(scalar);
		return a;
	}

	template <typename Scalar, class PyClass, typename boost::enable_if_c<Eigen::NumTraits<Scalar>::IsComplex, int>::type = 0>
	static void visit_reductions_noncomplex(PyClass&)
	{ /* do nothing*/
	}
	template <typename Scalar, class PyClass, typename boost::disable_if_c<Eigen::NumTraits<Scalar>::IsComplex, int>::type = 0>
	static void visit_reductions_noncomplex(PyClass& cl)
	{
		// must be wrapped since there are overloads:
		//   maxCoeff(), maxCoeff(IndexType*), maxCoeff(IndexType*,IndexType*)
		cl.def("maxCoeff", &MatrixBaseVisitor::maxCoeff0, "Maximum value over all elements.")
		        .def("minCoeff", &MatrixBaseVisitor::minCoeff0, "Minimum value over all elements.");
	}
	static Scalar maxCoeff0(const MatrixBaseT& m) { return m.array().maxCoeff(); }
	static Scalar minCoeff0(const MatrixBaseT& m) { return m.array().minCoeff(); }

	// we want to keep -0 (rather than replacing it by 0), but that does not work for complex numbers
	// hence two versions
	template <typename Scalar, typename boost::disable_if_c<Eigen::NumTraits<Scalar>::IsComplex, int>::type = 0>
	static bool prune_element(const Scalar& num, RealScalar absTol)
	{
		using namespace std;
		return abs(num) <= absTol || num != -0;
	}
	template <typename Scalar, typename boost::enable_if_c<Eigen::NumTraits<Scalar>::IsComplex, int>::type = 0>
	static bool prune_element(const Scalar& num, RealScalar absTol)
	{
		using namespace std;
		return abs(num) <= absTol;
	}

	static MatrixBaseT pruned(const MatrixBaseT& a, double absTol = 1e-6)
	{ // typename MatrixBaseT::Scalar absTol=1e-6){
		MatrixBaseT ret(MatrixBaseT::Zero(a.rows(), a.cols()));
		for (Index c = 0; c < a.cols(); c++) {
			for (Index r = 0; r < a.rows(); r++) {
				if (!prune_element(a(c, r), absTol)) ret(c, r) = a(c, r);
			}
		}
		return ret;
	};
};

template <typename VectorT> class VectorVisitor : public py::def_visitor<VectorVisitor<VectorT>> {
	friend class def_visitor_access;
	typedef typename VectorT::Scalar                                                      Scalar;
	typedef typename ::yade::math::RealOf<Scalar>                                         RealT;
	typedef Eigen::Matrix<Scalar, VectorT::RowsAtCompileTime, VectorT::RowsAtCompileTime> CompatMatrixT;
	typedef Eigen::Matrix<Scalar, 2, 1>                                                   CompatVec2;
	typedef Eigen::Matrix<Scalar, 3, 1>                                                   CompatVec3;
	typedef Eigen::Matrix<Scalar, 4, 1>                                                   CompatVec4;
	typedef Eigen::Matrix<Scalar, 6, 1>                                                   CompatVec6;
	typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1>                                      CompatVecX;
	enum { Dim = VectorT::RowsAtCompileTime };

public:
	template <class PyClass> void visit(PyClass& cl) const
	{
		MatrixBaseVisitor<VectorT>().visit(cl);
		cl.def_pickle(VectorPickle())
		        .def("__setitem__", &VectorVisitor::set_item)
		        .def("__getitem__", &VectorVisitor::get_item)
		        .def("__str__", &VectorVisitor::__str__)
		        .def("__repr__", &VectorVisitor::__str__)
		        .def("dot", &VectorVisitor::dot, py::arg("other"), "Dot product with *other*.")
		        .def("outer", &VectorVisitor::outer, py::arg("other"), "Outer product with *other*.")
		        .def("asDiagonal", &VectorVisitor::asDiagonal, "Return diagonal matrix with this vector on the diagonal.");

		visit_fixed_or_dynamic<VectorT, PyClass>(cl);

		visit_special_sizes<VectorT, PyClass>(cl);
	};

private:
	// for dynamic vectors
	template <typename VectorT2, class PyClass, typename boost::enable_if_c<VectorT2::RowsAtCompileTime == Eigen::Dynamic, int>::type = 0>
	static void visit_fixed_or_dynamic(PyClass& cl)
	{
		//std::cerr<<"VectorVisitor: dynamic vector for "<<name()<<std::endl;
		cl.def("__len__", &VectorVisitor::dyn__len__)
		        .def("resize", &VectorVisitor::resize)
		        .def("Unit", &VectorVisitor::dyn_Unit)
		        .staticmethod("Unit")
		        .def("Ones", &VectorVisitor::dyn_Ones)
		        .staticmethod("Ones")
		        .def("Zero", &VectorVisitor::dyn_Zero)
		        .staticmethod("Zero")
		        .def("Random",
		             &VectorVisitor::dyn_Random,
		             py::arg("len"),
		             "Return vector of given length with all elements set to values between 0 and 1 randomly.")
		        .staticmethod("Random");
	}
	// for fixed-size vectors
	template <typename VectorT2, class PyClass, typename boost::disable_if_c<VectorT2::RowsAtCompileTime == Eigen::Dynamic, int>::type = 0>
	static void visit_fixed_or_dynamic(PyClass& cl)
	{
		//std::cerr<<"VectorVisitor: fixed vector for "<<name()<<std::endl;
		cl.def("__len__", &VectorVisitor::__len__).staticmethod("__len__").def("Unit", &VectorVisitor::Unit).staticmethod("Unit");
	}

	// handle specific sizes of vectors separately

	// 2-vector
	template <typename VectorT2, class PyClass, typename boost::enable_if_c<VectorT2::RowsAtCompileTime == 2, int>::type = 0>
	static void visit_special_sizes(PyClass& cl)
	{
		cl.def("__init__", py::make_constructor(&VectorVisitor::fromStrV2, py::default_call_policies(), (py::arg("str1"), py::arg("str2"))))
		        .def(py::init<typename VectorT2::Scalar, typename VectorT2::Scalar>((py::arg("x"), py::arg("y"))))
		        .add_static_property("UnitX", &VectorVisitor::Vec2_UnitX)
		        .add_static_property("UnitY", &VectorVisitor::Vec2_UnitY);
	}
	static CompatVec2* fromStrV2(const std::string& s1, const std::string& s2)
	{
		CompatVec2* v2(new CompatVec2(conv::val<Scalar>(s1), conv::val<Scalar>(s2)));
		return v2;
	}
	static CompatVec2 Vec2_UnitX() { return CompatVec2::UnitX(); }
	static CompatVec2 Vec2_UnitY() { return CompatVec2::UnitY(); }

	// 3-vector
	template <typename VectorT2, class PyClass, typename boost::enable_if_c<VectorT2::RowsAtCompileTime == 3, int>::type = 0>
	static void visit_special_sizes(PyClass& cl)
	{
		cl.def("__init__",
		       py::make_constructor(&VectorVisitor::fromStrV3, py::default_call_policies(), (py::arg("str1"), py::arg("str2"), py::arg("str3"))))
		        .def(py::init<typename VectorT2::Scalar, typename VectorT2::Scalar, typename VectorT2::Scalar>(
		                (py::arg("x") = Scalar(0), py::arg("y") = Scalar(0), py::arg("z") = Scalar(0))))
		        .def("cross", &VectorVisitor::cross) // cross-product only meaningful for 3-sized vectors
		        .add_static_property("UnitX", &VectorVisitor::Vec3_UnitX)
		        .add_static_property("UnitY", &VectorVisitor::Vec3_UnitY)
		        .add_static_property("UnitZ", &VectorVisitor::Vec3_UnitZ)
		        // swizzles
		        .def("xy", &VectorVisitor::Vec3_xy)
		        .def("yx", &VectorVisitor::Vec3_yx)
		        .def("xz", &VectorVisitor::Vec3_xz)
		        .def("zx", &VectorVisitor::Vec3_zx)
		        .def("yz", &VectorVisitor::Vec3_yz)
		        .def("zy", &VectorVisitor::Vec3_zy);
	}
	static CompatVec3* fromStrV3(const std::string& s1, const std::string& s2, const std::string& s3)
	{
		CompatVec3* v3(new CompatVec3(conv::val<Scalar>(s1), conv::val<Scalar>(s2), conv::val<Scalar>(s3)));
		return v3;
	}
	static CompatVec3 cross(const CompatVec3& self, const CompatVec3& other) { return self.cross(other); }
	static CompatVec3 Vec3_UnitX() { return CompatVec3::UnitX(); }
	static CompatVec3 Vec3_UnitY() { return CompatVec3::UnitY(); }
	static CompatVec3 Vec3_UnitZ() { return CompatVec3::UnitZ(); }

	static CompatVec2 Vec3_xy(const CompatVec3& v) { return CompatVec2(v[0], v[1]); }
	static CompatVec2 Vec3_yx(const CompatVec3& v) { return CompatVec2(v[1], v[0]); }
	static CompatVec2 Vec3_xz(const CompatVec3& v) { return CompatVec2(v[0], v[2]); }
	static CompatVec2 Vec3_zx(const CompatVec3& v) { return CompatVec2(v[2], v[0]); }
	static CompatVec2 Vec3_yz(const CompatVec3& v) { return CompatVec2(v[1], v[2]); }
	static CompatVec2 Vec3_zy(const CompatVec3& v) { return CompatVec2(v[2], v[1]); }

	// 4-vector
	template <typename VectorT2, class PyClass, typename boost::enable_if_c<VectorT2::RowsAtCompileTime == 4, int>::type = 0>
	static void visit_special_sizes(PyClass& cl)
	{
		cl.def("__init__",
		       py::make_constructor(
		               &VectorVisitor::fromStrV4, py::default_call_policies(), (py::arg("str1"), py::arg("str2"), py::arg("str3"), py::arg("str4"))))
		        .def(py::init<typename VectorT2::Scalar, typename VectorT2::Scalar, typename VectorT2::Scalar, typename VectorT2::Scalar>(
		                (py::arg("v0"), py::arg("v1"), py::arg("v2"), py::arg("v3"))));
	}
	static CompatVec4* fromStrV4(const std::string& s1, const std::string& s2, const std::string& s3, const std::string& s4)
	{
		CompatVec4* v4(new CompatVec4(conv::val<Scalar>(s1), conv::val<Scalar>(s2), conv::val<Scalar>(s3), conv::val<Scalar>(s4)));
		return v4;
	}

	// 6-vector
	template <typename VectorT2, class PyClass, typename boost::enable_if_c<VectorT2::RowsAtCompileTime == 6, int>::type = 0>
	static void visit_special_sizes(PyClass& cl)
	{
		cl.def("__init__",
		       py::make_constructor(
		               &VectorVisitor::Vec6_fromElements,
		               py::default_call_policies(),
		               (py::arg("v0"), py::arg("v1"), py::arg("v2"), py::arg("v3"), py::arg("v4"), py::arg("v5"))))
		        .def("__init__",
		             py::make_constructor(&VectorVisitor::Vec6_fromHeadTail, py::default_call_policies(), (py::arg("head"), py::arg("tail"))))
		        .def("head", &VectorVisitor::Vec6_head)
		        .def("tail", &VectorVisitor::Vec6_tail);
	}
	// only used with dim==6
	static CompatVec6* Vec6_fromElements(const Scalar& v0, const Scalar& v1, const Scalar& v2, const Scalar& v3, const Scalar& v4, const Scalar& v5)
	{
		CompatVec6* v(new CompatVec6);
		(*v) << v0, v1, v2, v3, v4, v5;
		return v;
	}
	static CompatVec6* Vec6_fromHeadTail(const CompatVec3& head, const CompatVec3& tail)
	{
		CompatVec6* v(new CompatVec6);
		(*v) << head, tail;
		return v;
	}
	static CompatVec3 Vec6_head(const CompatVec6& v) { return v.template head<3>(); }
	static CompatVec3 Vec6_tail(const CompatVec6& v) { return v.template tail<3>(); }

	// ctor for dynamic vectors
	template <typename VectorT2, class PyClass, typename boost::enable_if_c<VectorT2::RowsAtCompileTime == Eigen::Dynamic, int>::type = 0>
	static void visit_special_sizes(PyClass& cl)
	{
		cl.def("__init__", py::make_constructor(&VecX_fromList, py::default_call_policies(), (py::arg("vv"))));
	}
	static CompatVecX* VecX_fromList(const std::vector<Scalar>& ii)
	{
		CompatVecX* v(new CompatVecX(ii.size()));
		for (size_t i = 0; i < ii.size(); i++)
			(*v)[i] = ii[i];
		return v;
	}


	static VectorT dyn_Ones(Index size) { return VectorT::Ones(size); }
	static VectorT dyn_Zero(Index size) { return VectorT::Zero(size); }
	static VectorT dyn_Random(Index size) { return VectorT::Random(size); }
	static VectorT Unit(Index ix)
	{
		IDX_CHECK(ix, (Index)Dim);
		return VectorT::Unit(ix);
	}
	static VectorT dyn_Unit(Index size, Index ix)
	{
		IDX_CHECK(ix, size);
		return VectorT::Unit(size, ix);
	}
	static bool  dyn() { return Dim == Eigen::Dynamic; }
	static Index __len__()
	{
		assert(!dyn());
		return Dim;
	}
	static Index         dyn__len__(const VectorT& self) { return self.size(); }
	static void          resize(VectorT& self, Index size) { self.resize(size); }
	static Scalar        dot(const VectorT& self, const VectorT& other) { return self.dot(other); }
	static CompatMatrixT outer(const VectorT& self, const VectorT& other) { return self * other.transpose(); }
	static CompatMatrixT asDiagonal(const VectorT& self) { return self.asDiagonal(); }
	static Scalar        get_item(const VectorT& self, Index ix)
	{
		IDX_CHECK(ix, dyn() ? (Index)self.size() : (Index)Dim);
		return self[ix];
	}
	static void set_item(VectorT& self, Index ix, Scalar value)
	{
		IDX_CHECK(ix, dyn() ? (Index)self.size() : (Index)Dim);
		self[ix] = value;
	}
	struct VectorPickle : py::pickle_suite {
		static py::tuple getinitargs(const VectorT& x)
		{
			// if this fails, add supported size to the switch below
			BOOST_STATIC_ASSERT(Dim == 2 || Dim == 3 || Dim == 4 || Dim == 6 || Dim == Eigen::Dynamic);
			switch ((Index)Dim) {
				case 2: return py::make_tuple(x[0], x[1]);
				case 3: return py::make_tuple(x[0], x[1], x[2]);
				case 4: return py::make_tuple(x[0], x[1], x[2], x[3]);
				case 6: return py::make_tuple(x[0], x[1], x[2], x[3], x[4], x[5]);
				default: return py::make_tuple(py::list(x));
			}
		};
	};

public:
	static std::string __str__(const py::object& obj)
	{
		std::ostringstream oss;
		const auto         self = py::extract<VectorT>(obj)();
		bool               list = (Dim == Eigen::Dynamic && self.size() > 0);
		oss << object_class_name(obj) << (list ? "([" : "(");
		Vector_data_stream(self, oss);
		oss << (list ? "])" : ")");
		return oss.str();
	};

	// not sure why this must be templated now?!
	template <typename VectorType> static void Vector_data_stream(const VectorType& self, std::ostringstream& oss, int pad = 0)
	{
		for (Index i = 0; i < self.size(); i++)
			oss << (i == 0 ? "" : (((i % 3) != 0 || pad > 0) ? "," : ", "))
			    << ::yade::minieigenHP::numToStringHP(self.row(i / self.cols())[i % self.cols()] /*,pad*/);
	}
};

template <typename MatrixT> class MatrixVisitor : public py::def_visitor<MatrixVisitor<MatrixT>> {
	friend class def_visitor_access;
	typedef typename MatrixT::Scalar                                      Scalar;
	typedef typename Eigen::Matrix<Scalar, MatrixT::RowsAtCompileTime, 1> CompatVectorT;
	typedef Eigen::Matrix<Scalar, 3, 3>                                   CompatMat3;
	typedef Eigen::Matrix<Scalar, 3, 1>                                   CompatVec3;
	typedef Eigen::Matrix<Scalar, 6, 6>                                   CompatMat6;
	typedef Eigen::Matrix<Scalar, 6, 1>                                   CompatVec6;
	typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>         CompatMatX;
	typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1>                      CompatVecX;
	enum { Dim = MatrixT::RowsAtCompileTime };

public:
	template <class PyClass> void visit(PyClass& cl) const
	{
		MatrixBaseVisitor<MatrixT>().visit(cl);
		cl.def_pickle(MatrixPickle())
		        .def("__init__", py::make_constructor(&MatrixVisitor::fromDiagonal, py::default_call_policies(), (py::arg("diag"))))

		        .def("determinant", &MatrixT::determinant, "Return matrix determinant.")
		        .def("trace", &MatrixT::trace, "Return sum of diagonal elements.")
		        .def("transpose", &MatrixVisitor::transpose, "Return transposed matrix.")
		        .def("diagonal", &MatrixVisitor::diagonal, "Return diagonal as vector.")
		        .def("row", &MatrixVisitor::row, py::arg("row"), "Return row as vector.")
		        .def("col", &MatrixVisitor::col, py::arg("col"), "Return column as vector.")
		        // matrix*matrix product
		        .def("__mul__", &MatrixVisitor::__mul__)
		        .def("__imul__", &MatrixVisitor::__imul__)
		        // matrix*vector product
		        .def("__mul__", &MatrixVisitor::__mul__vec)
		        .def("__rmul__", &MatrixVisitor::__mul__vec)
		        .def("__setitem__", &MatrixVisitor::set_row)
		        .def("__getitem__", &MatrixVisitor::get_row)
		        .def("__setitem__", &MatrixVisitor::set_item)
		        .def("__getitem__", &MatrixVisitor::get_item)
		        .def("__str__", &MatrixVisitor::__str__)
		        .def("__repr__", &MatrixVisitor::__str__);
		visit_if_float<Scalar, PyClass>(cl);
		visit_fixed_or_dynamic<MatrixT, PyClass>(cl);
		visit_special_sizes<MatrixT, PyClass>(cl);

		//std::cerr<<"MatrixVisitor: "<<name()<<std::endl;
	}

private:
	// for dynamic matrices
	template <typename MatrixT2, class PyClass, typename boost::enable_if_c<MatrixT2::RowsAtCompileTime == Eigen::Dynamic, int>::type = 0>
	static void visit_fixed_or_dynamic(PyClass& cl)
	{
		cl.def("__len__", &MatrixVisitor::dyn__len__)
		        .def("resize",
		             &MatrixVisitor::resize,
		             "Change size of the matrix, keep values of elements which exist in the new matrix",
		             (py::arg("rows"), py::arg("cols")))
		        .def("Ones",
		             &MatrixVisitor::dyn_Ones,
		             (py::arg("rows"), py::arg("cols")),
		             "Create matrix of given dimensions where all elements are set to 1.")
		        .staticmethod("Ones")
		        .def("Zero", &MatrixVisitor::dyn_Zero, (py::arg("rows"), py::arg("cols")), "Create zero matrix of given dimensions")
		        .staticmethod("Zero")
		        .def("Random",
		             &MatrixVisitor::dyn_Random,
		             (py::arg("rows"), py::arg("cols")),
		             "Create matrix with given dimensions where all elements are set to number between 0 and 1 (uniformly-distributed).")
		        .staticmethod("Random")
		        .def("Identity", &MatrixVisitor::dyn_Identity, (py::arg("rank")), "Create identity matrix with given rank (square).")
		        .staticmethod("Identity");
	}
	// for fixed-size matrices
	template <typename MatrixT2, class PyClass, typename boost::disable_if_c<MatrixT2::RowsAtCompileTime == Eigen::Dynamic, int>::type = 0>
	static void visit_fixed_or_dynamic(PyClass& cl)
	{
		cl.def("__len__", &MatrixVisitor::__len__).staticmethod("__len__");
	}

	template <typename Scalar, class PyClass, typename boost::enable_if<boost::is_integral<Scalar>, int>::type = 0> static void visit_if_float(PyClass&)
	{ /* do nothing */
	}
	template <typename Scalar, class PyClass, typename boost::disable_if<boost::is_integral<Scalar>, int>::type = 0> static void visit_if_float(PyClass& cl)
	{
		cl
		        // matrix-matrix division?!
		        //.def("__div__",&MatrixBaseVisitor::__div__).def("__idiv__",&MatrixBaseVisitor::__idiv__)
		        .def("inverse", &MatrixVisitor::inverse, "Return inverted matrix.");
		// decompositions are only meaningful on non-complex numbers
		visit_if_decompositions_meaningful<Scalar, PyClass>(cl);
	}
	// for complex numbers, do nothing
	template <typename Scalar, class PyClass, typename boost::enable_if_c<Eigen::NumTraits<Scalar>::IsComplex, int>::type = 0>
	static void visit_if_decompositions_meaningful(PyClass&)
	{ /* do nothing */
	}
	// for non-complex numbers, define decompositions
	template <typename Scalar, class PyClass, typename boost::disable_if_c<Eigen::NumTraits<Scalar>::IsComplex, int>::type = 0>
	static void visit_if_decompositions_meaningful(PyClass& cl)
	{
		cl.def("jacobiSVD", &MatrixVisitor::jacobiSVD, "Compute SVD decomposition of square matrix, retuns (U,S,V) such that self=U*S*V.transpose()")
		        .def("svd", &MatrixVisitor::jacobiSVD, "Alias for :obj:`jacobiSVD`.")
		        .def("computeUnitaryPositive",
		             &MatrixVisitor::computeUnitaryPositive,
		             "Compute polar decomposition (unitary matrix U and positive semi-definite symmetric matrix P such that self=U*P).")
		        .def("polarDecomposition", &MatrixVisitor::computeUnitaryPositive, "Alias for :obj:`computeUnitaryPositive`.")
		        .def("selfAdjointEigenDecomposition",
		             &MatrixVisitor::selfAdjointEigenDecomposition,
		             "Compute eigen (spectral) decomposition of symmetric matrix, returns (eigVecs,eigVals). eigVecs is orthogonal Matrix3 with "
		             "columns ar normalized eigenvectors, eigVals is Vector3 with corresponding eigenvalues. "
		             "self=eigVecs*diag(eigVals)*eigVecs.transpose().")
		        .def("spectralDecomposition", &MatrixVisitor::selfAdjointEigenDecomposition, "Alias for :obj:`selfAdjointEigenDecomposition`.");
	}

	// handle specific matrix sizes
	// 3x3
	template <typename MatT2, class PyClass, typename boost::enable_if_c<MatT2::RowsAtCompileTime == 3, int>::type = 0>
	static void visit_special_sizes(PyClass& cl)
	{
		cl.def("__init__",
		       py::make_constructor(
		               &MatrixVisitor::Mat3_fromElements,
		               py::default_call_policies(),
		               (py::arg("m00"),
		                py::arg("m01"),
		                py::arg("m02"),
		                py::arg("m10"),
		                py::arg("m11"),
		                py::arg("m12"),
		                py::arg("m20"),
		                py::arg("m21"),
		                py::arg("m22"))))
		        .def("__init__",
		             py::make_constructor(
		                     &MatrixVisitor::Mat3_fromElementStrings,
		                     py::default_call_policies(),
		                     (py::arg("m00"),
		                      py::arg("m01"),
		                      py::arg("m02"),
		                      py::arg("m10"),
		                      py::arg("m11"),
		                      py::arg("m12"),
		                      py::arg("m20"),
		                      py::arg("m21"),
		                      py::arg("m22"))))
		        .def("__init__",
		             py::make_constructor(
		                     &MatrixVisitor::Mat3_fromRows,
		                     py::default_call_policies(),
		                     (py::arg("r0"), py::arg("r1"), py::arg("r2"), py::arg("cols") = false)));
	}
	static CompatMat3* Mat3_fromElements(
	        const Scalar& m00,
	        const Scalar& m01,
	        const Scalar& m02,
	        const Scalar& m10,
	        const Scalar& m11,
	        const Scalar& m12,
	        const Scalar& m20,
	        const Scalar& m21,
	        const Scalar& m22)
	{
		CompatMat3* m(new CompatMat3);
		(*m) << m00, m01, m02, m10, m11, m12, m20, m21, m22;
		return m;
	}
	static CompatMat3* Mat3_fromElementStrings(
	        const std::string& m00,
	        const std::string& m01,
	        const std::string& m02,
	        const std::string& m10,
	        const std::string& m11,
	        const std::string& m12,
	        const std::string& m20,
	        const std::string& m21,
	        const std::string& m22)
	{
		CompatMat3* m(new CompatMat3);
		(*m) << conv::val<Scalar>(m00), conv::val<Scalar>(m01), conv::val<Scalar>(m02), conv::val<Scalar>(m10), conv::val<Scalar>(m11),
		        conv::val<Scalar>(m12), conv::val<Scalar>(m20), conv::val<Scalar>(m21), conv::val<Scalar>(m22);
		return m;
	}
	static CompatMat3* Mat3_fromRows(const CompatVec3& l0, const CompatVec3& l1, const CompatVec3& l2, bool cols = false)
	{
		CompatMat3* m(new CompatMat3);
		if (cols) {
			m->col(0) = l0;
			m->col(1) = l1;
			m->col(2) = l2;
		} else {
			m->row(0) = l0;
			m->row(1) = l1;
			m->row(2) = l2;
		}
		return m;
	}

	// 6x6
	template <typename MatT2, class PyClass, typename boost::enable_if_c<MatT2::RowsAtCompileTime == 6, int>::type = 0>
	static void visit_special_sizes(PyClass& cl)
	{
		cl.def("__init__",
		       py::make_constructor(
		               &MatrixVisitor::Mat6_fromBlocks, py::default_call_policies(), (py::arg("ul"), py::arg("ur"), py::arg("ll"), py::arg("lr"))))
		        .def("__init__",
		             py::make_constructor(
		                     &MatrixVisitor::Mat6_fromRows,
		                     py::default_call_policies(),
		                     (py::arg("l0"), py::arg("l1"), py::arg("l2"), py::arg("l3"), py::arg("l4"), py::arg("l5"), py::arg("cols") = false)))
		        /* 3x3 blocks */
		        .def("ul", &MatrixVisitor::Mat6_ul, "Return upper-left 3x3 block")
		        .def("ur", &MatrixVisitor::Mat6_ur, "Return upper-right 3x3 block")
		        .def("ll", &MatrixVisitor::Mat6_ll, "Return lower-left 3x3 block")
		        .def("lr", &MatrixVisitor::Mat6_lr, "Return lower-right 3x3 block");
	}
	static CompatMat6* Mat6_fromBlocks(const CompatMat3& ul, const CompatMat3& ur, const CompatMat3& ll, const CompatMat3& lr)
	{
		CompatMat6* m(new CompatMat6);
		(*m) << ul, ur, ll, lr;
		return m;
	}
	static CompatMat6* Mat6_fromRows(
	        const CompatVec6& l0,
	        const CompatVec6& l1,
	        const CompatVec6& l2,
	        const CompatVec6& l3,
	        const CompatVec6& l4,
	        const CompatVec6& l5,
	        bool              cols = false)
	{
		CompatMat6* m(new CompatMat6);
		if (cols) {
			m->col(0) = l0;
			m->col(1) = l1;
			m->col(2) = l2;
			m->col(3) = l3;
			m->col(4) = l4;
			m->col(5) = l5;
		} else {
			m->row(0) = l0;
			m->row(1) = l1;
			m->row(2) = l2;
			m->row(3) = l3;
			m->row(4) = l4;
			m->row(5) = l5;
		}
		return m;
	}
	// see http://stackoverflow.com/questions/20870648/error-in-seemingly-correct-template-code-whats-wrong#20870661
	// on why the template keyword is needed after the m.
	static CompatMat3 Mat6_ul(const CompatMat6& m) { return m.template topLeftCorner<3, 3>(); }
	static CompatMat3 Mat6_ur(const CompatMat6& m) { return m.template topRightCorner<3, 3>(); }
	static CompatMat3 Mat6_ll(const CompatMat6& m) { return m.template bottomLeftCorner<3, 3>(); }
	static CompatMat3 Mat6_lr(const CompatMat6& m) { return m.template bottomRightCorner<3, 3>(); }

	// XxX
	template <typename MatT2, class PyClass, typename boost::enable_if_c<MatT2::RowsAtCompileTime == Eigen::Dynamic, int>::type = 0>
	static void visit_special_sizes(PyClass& cl)
	{
		cl.def("__init__",
		       py::make_constructor(
		               &MatrixVisitor::MatX_fromRows,
		               py::default_call_policies(),
		               (py::arg("r0")   = CompatVecX(),
		                py::arg("r1")   = CompatVecX(),
		                py::arg("r2")   = CompatVecX(),
		                py::arg("r3")   = CompatVecX(),
		                py::arg("r4")   = CompatVecX(),
		                py::arg("r5")   = CompatVecX(),
		                py::arg("r6")   = CompatVecX(),
		                py::arg("r7")   = CompatVecX(),
		                py::arg("r8")   = CompatVecX(),
		                py::arg("r9")   = CompatVecX(),
		                py::arg("cols") = false)))
		        .def("__init__",
		             py::make_constructor(&MatrixVisitor::MatX_fromRowSeq, py::default_call_policies(), (py::arg("rows"), py::arg("cols") = false)));
	}

	static CompatMatX* MatX_fromRows(
	        const CompatVecX& r0,
	        const CompatVecX& r1,
	        const CompatVecX& r2,
	        const CompatVecX& r3,
	        const CompatVecX& r4,
	        const CompatVecX& r5,
	        const CompatVecX& r6,
	        const CompatVecX& r7,
	        const CompatVecX& r8,
	        const CompatVecX& r9,
	        bool              setCols)
	{
		/* check vector dimensions */ CompatVecX rr[] = { r0, r1, r2, r3, r4, r5, r6, r7, r8, r9 };
		int                                      cols = -1, rows = -1;
		for (int i = 0; i < 10; i++) {
			if (rows < 0 && rr[i].size() == 0) rows = i;
			if (rows >= 0 && rr[i].size() > 0)
				throw std::invalid_argument("Matrix6r: non-empty rows not allowed after first empty row, which marks end of the matrix.");
		}
		cols = (rows > 0 ? rr[0].size() : 0);
		for (int i = 1; i < rows; i++)
			if (rr[i].size() != cols)
				throw std::invalid_argument(("Matrix6: all non-empty rows must have the same length (0th row has "
				                             + boost::lexical_cast<std::string>(rr[0].size()) + " items, " + boost::lexical_cast<std::string>(i)
				                             + "th row has " + boost::lexical_cast<std::string>(rr[i].size()) + " items)")
				                                    .c_str());
		CompatMatX* m;
		m = setCols ? new CompatMatX(cols, rows) : new CompatMatX(rows, cols);
		for (int i = 0; i < rows; i++) {
			if (setCols) m->col(i) = rr[i];
			else
				m->row(i) = rr[i];
		}
		return m;
	}
	static CompatMatX* MatX_fromRowSeq(const std::vector<CompatVecX>& rr, bool setCols)
	{
		int rows = rr.size(), cols = rr.size() > 0 ? rr[0].size() : 0;
		for (int i = 1; i < rows; i++)
			if (rr[i].size() != cols) throw std::invalid_argument(("MatrixX: all rows must have the same length."));
		CompatMatX* m;
		m = setCols ? new CompatMatX(cols, rows) : new CompatMatX(rows, cols);
		for (int i = 0; i < rows; i++) {
			if (setCols) m->col(i) = rr[i];
			else
				m->row(i) = rr[i];
		}
		return m;
	};


	static MatrixT                 dyn_Ones(Index rows, Index cols) { return MatrixT::Ones(rows, cols); }
	static MatrixT                 dyn_Zero(Index rows, Index cols) { return MatrixT::Zero(rows, cols); }
	static MatrixT                 dyn_Random(Index rows, Index cols) { return MatrixT::Random(rows, cols); }
	static MatrixT                 dyn_Identity(Index rows, Index cols) { return MatrixT::Identity(rows, cols); }
	static typename MatrixT::Index dyn__len__(MatrixT& a) { return a.rows(); }
	static typename MatrixT::Index __len__() { return MatrixT::RowsAtCompileTime; }
	static MatrixT                 Identity() { return MatrixT::Identity(); }
	static MatrixT                 transpose(const MatrixT& m) { return m.transpose(); }
	static CompatVectorT           diagonal(const MatrixT& m) { return m.diagonal(); }
	static MatrixT*                fromDiagonal(const CompatVectorT& d)
	{
		MatrixT* m(new MatrixT);
		*m = d.asDiagonal();
		return m;
	}
	static void          resize(MatrixT& self, Index rows, Index cols) { self.resize(rows, cols); }
	static CompatVectorT get_row(const MatrixT& a, Index ix)
	{
		IDX_CHECK(ix, a.rows());
		return a.row(ix);
	}
	static void set_row(MatrixT& a, Index ix, const CompatVectorT& r)
	{
		IDX_CHECK(ix, a.rows());
		a.row(ix) = r;
	}
	static Scalar get_item(const MatrixT& a, py::tuple _idx)
	{
		Index idx[2];
		Index mx[2] = { a.rows(), a.cols() };
		IDX2_CHECKED_TUPLE_INTS(_idx, mx, idx);
		return a(idx[0], idx[1]);
	}
	static void set_item(MatrixT& a, py::tuple _idx, const Scalar& value)
	{
		Index idx[2];
		Index mx[2] = { a.rows(), a.cols() };
		IDX2_CHECKED_TUPLE_INTS(_idx, mx, idx);
		a(idx[0], idx[1]) = value;
	}

	static MatrixT __imul__(MatrixT& a, const MatrixT& b)
	{
		a *= b;
		return a;
	};
	static MatrixT       __mul__(const MatrixT& a, const MatrixT& b) { return a * b; }
	static CompatVectorT __mul__vec(const MatrixT& m, const CompatVectorT& v) { return m * v; }
	// float matrices only
	static MatrixT inverse(const MatrixT& m) { return m.inverse(); }
	static MatrixT __div__(const MatrixT& a, const MatrixT& b) { return a / b; }
	// static void __idiv__(MatrixT& a, const MatrixT& b){ a/=b; };
	static CompatVectorT row(const MatrixT& m, Index ix)
	{
		IDX_CHECK(ix, m.rows());
		return m.row(ix);
	}
	static CompatVectorT col(const MatrixT& m, Index ix)
	{
		IDX_CHECK(ix, m.cols());
		return m.col(ix);
	}

	static void ensureSquare(const MatrixT& m)
	{
		if (m.rows() != m.cols()) throw std::runtime_error("Matrix is not square.");
	}
	static py::tuple jacobiSVD(const MatrixT& in)
	{
		ensureSquare(in);
		Eigen::JacobiSVD<MatrixT> svd(in, Eigen::ComputeThinU | Eigen::ComputeThinV);
		return py::make_tuple(svd.matrixU(), svd.matrixV(), MatrixT(svd.singularValues().asDiagonal()));
	};
	// polar decomposition
	static py::tuple computeUnitaryPositive(const MatrixT& in)
	{
		ensureSquare(in);
		Eigen::JacobiSVD<MatrixT> svd(in, Eigen::ComputeThinU | Eigen::ComputeThinV);
		const MatrixT&            u = svd.matrixU();
		const MatrixT&            v = svd.matrixV();
		MatrixT                   s = svd.singularValues().asDiagonal();
		return py::make_tuple(u * v.transpose(), v * s * v.transpose());
	}
	// eigen decomposition
	static py::tuple selfAdjointEigenDecomposition(const MatrixT& in)
	{
		ensureSquare(in);
		Eigen::SelfAdjointEigenSolver<MatrixT> a(in);
		return py::make_tuple(a.eigenvectors(), a.eigenvalues());
	}
	static bool        dyn() { return Dim == Eigen::Dynamic; }
	static std::string __str__(const py::object& obj)
	{
		std::ostringstream oss;
		const auto         m = py::extract<MatrixT>(obj)();
		oss << object_class_name(obj) << "(";
		bool wrap = ((dyn() && m.rows() > 1) || (!dyn() && m.rows() > 3));
		// non-wrapping fixed-size: flat list of numbers, not rows as tuples (Matrix3)
		if (!dyn() && !wrap) {
			VectorVisitor<CompatVectorT>::template Vector_data_stream<MatrixT>(m, oss, /*pad=*/0);
		} else {
			if (wrap) oss << "\n";
			for (Index r = 0; r < m.rows(); r++) {
				oss << (wrap ? "\t" : "") << "(";
				VectorVisitor<CompatVectorT>::template Vector_data_stream<CompatVectorT>(m.row(r), oss, /*pad=*/(wrap ? 7 : 0));
				oss << ")" << (r < m.rows() - 1 ? "," : "") << (wrap ? "\n" : "");
			}
		}
		oss << ")";
		return oss.str();
	}
	struct MatrixPickle : py::pickle_suite {
		static py::tuple getinitargs(const MatrixT& x)
		{
			// if this fails, add supported size to the switch below
			BOOST_STATIC_ASSERT(Dim == 2 || Dim == 3 || Dim == 6 || Dim == Eigen::Dynamic);
			switch ((Index)Dim) {
				case 2: return py::make_tuple(x(0, 0), x(0, 1), x(1, 0), x(1, 1));
				case 3: return py::make_tuple(x(0, 0), x(0, 1), x(0, 2), x(1, 0), x(1, 1), x(1, 2), x(2, 0), x(2, 1), x(2, 2));
				case 6: return py::make_tuple(x.row(0), x.row(1), x.row(2), x.row(3), x.row(4), x.row(5));
				// should return list of rows, which are VectorX
				default: return py::make_tuple(py::list(x));
			}
		};
	};
};


template <typename Box> class AabbVisitor : public py::def_visitor<AabbVisitor<Box>> {
	friend class def_visitor_access;
	typedef typename Box::VectorType VectorType;
	typedef typename Box::Scalar     Scalar;

public:
	template <class PyClass> void visit(PyClass& cl) const
	{
		cl.def(py::init<Box>(py::arg("other")))
		        .def(py::init<VectorType, VectorType>((py::arg("min"), py::arg("max"))))
		        .def_pickle(BoxPickle())
		        .def("volume", &Box::volume)
		        .def("empty", &Box::isEmpty)
		        .def("center", &AabbVisitor::center)
		        .def("sizes", &AabbVisitor::sizes)
		        .def("contains", &AabbVisitor::containsPt)
		        .def("contains", &AabbVisitor::containsBox)
		        // for the "in" operator
		        .def("__contains__", &AabbVisitor::containsPt)
		        .def("__contains__", &AabbVisitor::containsBox)
		        .def("extend", &AabbVisitor::extendPt)
		        .def("extend", &AabbVisitor::extendBox)
		        .def("clamp", &AabbVisitor::clamp)
		        // return new objects
		        .def("intersection", &Box::intersection)
		        .def("merged", &Box::merged)
		        // those return internal references, which is what we want
		        .add_property("min", py::make_function(&AabbVisitor::min, py::return_internal_reference<>()))
		        .add_property("max", py::make_function(&AabbVisitor::max, py::return_internal_reference<>()))
		        .def("__len__", &AabbVisitor::len)
		        .staticmethod("__len__")
		        .def("__setitem__", &AabbVisitor::set_item)
		        .def("__getitem__", &AabbVisitor::get_item)
		        .def("__setitem__", &AabbVisitor::set_minmax)
		        .def("__getitem__", &AabbVisitor::get_minmax)
		        .def("__str__", &AabbVisitor::__str__)
		        .def("__repr__", &AabbVisitor::__str__);
	};

private:
	static bool        containsPt(const Box& self, const VectorType& pt) { return self.contains(pt); }
	static bool        containsBox(const Box& self, const Box& other) { return self.contains(other); }
	static void        extendPt(Box& self, const VectorType& pt) { self.extend(pt); }
	static void        extendBox(Box& self, const Box& other) { self.extend(other); }
	static void        clamp(Box& self, const Box& other) { self.clamp(other); }
	static VectorType& min(Box& self) { return self.min(); }
	static VectorType& max(Box& self) { return self.max(); }
	static VectorType  center(const Box& self) { return self.center(); }
	static VectorType  sizes(const Box& self) { return self.sizes(); }
	struct BoxPickle : py::pickle_suite {
		static py::tuple getinitargs(const Box& x) { return py::make_tuple(x.min(), x.max()); }
	};
	static Index len() { return 2; }
	// getters and setters
	static Scalar get_item(const Box& self, py::tuple _idx)
	{
		Index idx[2];
		Index mx[2] = { 2, Box::AmbientDimAtCompileTime };
		IDX2_CHECKED_TUPLE_INTS(_idx, mx, idx);
		if (idx[0] == 0) return self.min()[idx[1]];
		return self.max()[idx[1]];
	}
	static void set_item(Box& self, py::tuple _idx, Scalar value)
	{
		Index idx[2];
		Index mx[2] = { 2, Box::AmbientDimAtCompileTime };
		IDX2_CHECKED_TUPLE_INTS(_idx, mx, idx);
		if (idx[0] == 0) self.min()[idx[1]] = value;
		else
			self.max()[idx[1]] = value;
	}
	static VectorType get_minmax(const Box& self, Index idx)
	{
		IDX_CHECK(idx, 2);
		if (idx == 0) return self.min();
		return self.max();
	}
	static void set_minmax(Box& self, Index idx, const VectorType& value)
	{
		IDX_CHECK(idx, 2);
		if (idx == 0) self.min() = value;
		else
			self.max() = value;
	}
	static std::string __str__(const py::object& obj)
	{
		const Box          self = py::extract<Box>(obj)();
		std::ostringstream oss;
		oss << object_class_name(obj) << "((";
		VectorVisitor<VectorType>::template Vector_data_stream<VectorType>(self.min(), oss);
		oss << "), (";
		VectorVisitor<VectorType>::template Vector_data_stream<VectorType>(self.max(), oss);
		oss << "))";
		return oss.str();
	}
};

template <typename QuaternionT, int Level = ::yade::math::levelOfRealHP<typename QuaternionT::Scalar>>
class QuaternionVisitor : public py::def_visitor<QuaternionVisitor<QuaternionT>> {
	typedef typename QuaternionT::Scalar             Scalar;
	typedef Eigen::Matrix<Scalar, 3, 1>              CompatVec3;
	typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> CompatVecX;
	typedef Eigen::Matrix<Scalar, 3, 3>              CompatMat3;
	typedef Eigen::AngleAxis<Scalar>                 AngleAxisT;

public:
	template <class PyClass> void visit(PyClass& cl) const
	{
		cl.def("__init__", py::make_constructor(&QuaternionVisitor::fromAxisAngleMpf, py::default_call_policies(), (py::arg("axis"), py::arg("angle"))))
		        .def("__init__",
		             py::make_constructor(&QuaternionVisitor::fromAxisAngle, py::default_call_policies(), (py::arg("axis"), py::arg("angle"))))
		        .def("__init__",
		             py::make_constructor(&QuaternionVisitor::fromAngleAxisMpf, py::default_call_policies(), (py::arg("angle"), py::arg("axis"))))
		        .def("__init__",
		             py::make_constructor(&QuaternionVisitor::fromAngleAxis, py::default_call_policies(), (py::arg("angle"), py::arg("axis"))))
		        .def("__init__",
		             py::make_constructor(&QuaternionVisitor::fromTupleStr, py::default_call_policies(), (py::arg("axis"), py::arg("angle"))))
		        .def("__init__", py::make_constructor(&QuaternionVisitor::fromTuple, py::default_call_policies(), (py::arg("tuple"))))
		        .def("__init__", py::make_constructor(&QuaternionVisitor::fromTwoVectors, py::default_call_policies(), (py::arg("u"), py::arg("v"))))
		        .def("__init__",
		             py::make_constructor(
		                     &QuaternionVisitor::fromStr,
		                     py::default_call_policies(),
		                     (py::arg("str1"), py::arg("str2"), py::arg("str3"), py::arg("str4"))))
		        .def(py::init<Scalar, Scalar, Scalar, Scalar>(
		                (py::arg("w"), py::arg("x"), py::arg("y"), py::arg("z")),
		                "Initialize from coefficients.\n\n.. note:: The order of coefficients is *w*, *x*, *y*, *z*. The [] operator numbers them "
		                "differently, 0...4 for *x* *y* *z* *w*!"))
		        .def(py::init<CompatMat3>((py::arg("rotMatrix")))) //,"Initialize from given rotation matrix.")
		        .def(py::init<QuaternionT>((py::arg("other"))))
		        .def_pickle(QuaternionPickle())
		        // properties
		        .add_static_property("Identity", &QuaternionVisitor::Identity)
		        // methods
		        .def("setFromTwoVectors", &QuaternionVisitor::setFromTwoVectors, ((py::arg("u"), py::arg("v"))))
		        .def("angularDistance", &QuaternionVisitor::angularDistance)
		        .def("conjugate", &QuaternionT::conjugate)
		        .def("toAxisAngle", &QuaternionVisitor::toAxisAngle)
		        .def("toAngleAxis", &QuaternionVisitor::toAngleAxis)
		        .def("toRotationMatrix", &QuaternionT::toRotationMatrix)
		        .def("toRotationVector", &QuaternionVisitor::toRotationVector)
		        .def("Rotate", &QuaternionVisitor::Rotate, ((py::arg("v"))))
		        .def("inverse", &QuaternionT::inverse)
		        .def("norm", &QuaternionT::norm)
		        .def("normalize", &QuaternionT::normalize)
		        .def("normalized", &QuaternionT::normalized)
		        .def("slerp", &QuaternionVisitor::slerp, (py::arg("t"), py::arg("other")))
		        // .def("random",&QuaternionVisitor::random,"Assign random orientation to the quaternion.")
		        // operators
		        .def(py::self * py::self)
		        .def(py::self *= py::self)
		        .def(py::self * py::other<CompatVec3>())
		        .def("__eq__", &QuaternionVisitor::__eq__)
		        .def("__ne__", &QuaternionVisitor::__ne__)
		        .def("__sub__", &QuaternionVisitor::__sub__)
		        // specials
		        .def("__abs__", &QuaternionT::norm)
		        .def("__len__", &QuaternionVisitor::__len__)
		        .staticmethod("__len__")
		        .def("__setitem__", &QuaternionVisitor::__setitem__)
		        .def("__getitem__", &QuaternionVisitor::__getitem__)
		        .def("__str__", &QuaternionVisitor::__str__)
		        .def("__repr__", &QuaternionVisitor::__str__);

		if (std::numeric_limits<Scalar>::digits10 >= 18) {
			cl.def(py::init<double, double, double, double>(
			        (py::arg("w"), py::arg("x"), py::arg("y"), py::arg("z")),
			        "Initialize from coefficients.\n\n.. note:: The order of coefficients is *w*, *x*, *y*, *z*. The [] operator numbers them "
			        "differently, 0...4 for *x* *y* *z* *w*!"));
		}
	}

private:
	static QuaternionT* fromAxisAngleMpf(const CompatVec3& axis, const py::object& angle) // mpmath support
	{
		QuaternionT* ret = new QuaternionT(AngleAxisT(::yade::math::fromStringRealHP<Scalar>(py::extract<std::string>(angle.attr("__str__")())), axis));
		ret->normalize();
		return ret;
	}
	static QuaternionT* fromAxisAngle(const CompatVec3& axis, const Scalar& angle)
	{
		QuaternionT* ret = new QuaternionT(AngleAxisT(angle, axis));
		ret->normalize();
		return ret;
	}
	static QuaternionT* fromAngleAxisMpf(const py::object& angle, const CompatVec3& axis) // mpmath support
	{
		QuaternionT* ret = new QuaternionT(AngleAxisT(::yade::math::fromStringRealHP<Scalar>(py::extract<std::string>(angle.attr("__str__")())), axis));
		ret->normalize();
		return ret;
	}
	static QuaternionT* fromAngleAxis(const Scalar& angle, const CompatVec3& axis)
	{
		QuaternionT* ret = new QuaternionT(AngleAxisT(angle, axis));
		ret->normalize();
		return ret;
	}
	static Scalar       S(const py::object& obj) { return ::yade::math::fromStringRealHP<Scalar>(py::extract<std::string>(obj.attr("__str__")())); }
	static Scalar       S(const std::string& s) { return ::yade::math::fromStringRealHP<Scalar>(s); }
	static QuaternionT* fromStr(const std::string& s1, const std::string& s2, const std::string& s3, const std::string& s4)
	{
		QuaternionT* q(new QuaternionT(S(s1), S(s2), S(s3), S(s4)));
		return q;
	}
	static QuaternionT* fromTupleStr(const py::tuple& axis, const std::string& angle)
	{
		return fromAxisAngle(CompatVec3(S(axis[0]), S(axis[1]), S(axis[2])), S(angle));
	}
	static QuaternionT* fromTuple(const py::tuple& tuple)
	{
		std::string s0 = py::extract<std::string>(tuple[0].attr("__class__").attr("__name__"));
		std::string s1 = py::extract<std::string>(tuple[1].attr("__class__").attr("__name__"));
		if (s0 == "Vector3") {
			if (s1 == "Real") { // Vector3,Real
				return fromAxisAngle(py::extract<CompatVec3>(tuple[0]), py::extract<Scalar>(tuple[1]));
			} else { // Vector3,mpmath.mpf
				return fromAxisAngleMpf(py::extract<CompatVec3>(tuple[0]), tuple[1]);
			}
		} else {
			if (s0 == "Real") { // Real, Vector3
				return fromAxisAngle(py::extract<CompatVec3>(tuple[1]), py::extract<Scalar>(tuple[0]));
			} else { // mpmath.mpf, Vector3
				return fromAxisAngleMpf(py::extract<CompatVec3>(tuple[1]), tuple[0]);
			}
		}
	}
	static QuaternionT* fromTwoVectors(const CompatVec3& u, const CompatVec3& v)
	{
		QuaternionT* q(new QuaternionT);
		q->setFromTwoVectors(u, v);
		return q;
	}

	// those must be wrapped since "other" is declared as QuaternionBase<OtherDerived>; the type is then not inferred when using .def
	static QuaternionT slerp(const QuaternionT& self, const Scalar& t, const QuaternionT& other) { return self.slerp(t, other); }
	static Scalar      angularDistance(const QuaternionT& self, const QuaternionT& other) { return self.angularDistance(other); }

	struct QuaternionPickle : py::pickle_suite {
		static py::tuple getinitargs(const QuaternionT& x) { return py::make_tuple(x.w(), x.x(), x.y(), x.z()); }
	};
	static QuaternionT       Identity() { return QuaternionT::Identity(); }
	static Vector3rHP<Level> Rotate(const QuaternionT& self, const Vector3rHP<Level>& u) { return self * u; }
	static py::tuple         toAxisAngle(const QuaternionT& self)
	{
		AngleAxisT aa(self);
		return py::make_tuple(aa.axis(), aa.angle());
	}
	static py::tuple toAngleAxis(const QuaternionT& self)
	{
		AngleAxisT aa(self);
		return py::make_tuple(aa.angle(), aa.axis());
	}
	static CompatVec3 toRotationVector(const QuaternionT& self)
	{
		AngleAxisT aa(self);
		return aa.angle() * aa.axis();
	}
	static void setFromTwoVectors(QuaternionT& self, const CompatVec3& u, const CompatVec3& v) { self.setFromTwoVectors(u, v); /*return self;*/ }

	static bool       __eq__(const QuaternionT& u, const QuaternionT& v) { return u.x() == v.x() && u.y() == v.y() && u.z() == v.z() && u.w() == v.w(); }
	static bool       __ne__(const QuaternionT& u, const QuaternionT& v) { return !__eq__(u, v); }
	static CompatVecX __sub__(const QuaternionT& a, const QuaternionT& b)
	{
		CompatVecX r(4);
		r << a.w() - b.w(), a.x() - b.x(), a.y() - b.y(), a.z() - b.z();
		return r;
	}

	static Scalar __getitem__(const QuaternionT& self, Index idx)
	{
		IDX_CHECK(idx, 4);
		if (idx == 0) return self.x();
		if (idx == 1) return self.y();
		if (idx == 2) return self.z();
		return self.w();
	}
	static void __setitem__(QuaternionT& self, Index idx, Scalar value)
	{
		IDX_CHECK(idx, 4);
		if (idx == 0) self.x() = value;
		else if (idx == 1)
			self.y() = value;
		else if (idx == 2)
			self.z() = value;
		else if (idx == 3)
			self.w() = value;
	}
	static std::string __str__(const py::object& obj)
	{
		const QuaternionT self = py::extract<QuaternionT>(obj)();
		AngleAxisT        aa(self);
		return std::string(object_class_name(obj) + "((") + ::yade::minieigenHP::numToStringHP(aa.axis()[0]) + ","
		        + ::yade::minieigenHP::numToStringHP(aa.axis()[1]) + "," + ::yade::minieigenHP::numToStringHP(aa.axis()[2]) + "),"
		        + ::yade::minieigenHP::numToStringHP(aa.angle()) + ")";
	}
	static Index __len__() { return 4; }
};
